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We analyze the scaling behavior of the two smallest Lyapunov exponents for electrons propagating 
on two-dimensional lattices with energies within a very narrow interval around the chiral critical 
point at E = in the presence of a perpendicular random magnetic flux. By a numerical analysis of 
the energy and size dependence we confirm that the two smallest Lyapunov exponents are functions 
of a single parameter. The latter is given by InL/ ln^(_E), which is the ratio of the logarithm of the 
system width L to the logarithm of the correlation length £(E). Close to the chiral critical point and 
energy \E\ <g Eo, we find a logarithmically divergent energy dependence hi£(E) oc | ln(_Eo/|i?|)| 1/ ' 2 , 
where Eo is a characteristic energy scale. Our data are in agreement with the theoretical prediction 
of M. Fabrizio and C. Castelliani [Nucl. Phys. B 583, 542 (2000)] and resolve an inconsistency of 
previous numerical work. 
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I. INTRODUCTION 

The numerical determination of transport parameters 
for electrons propagating in disordered two-dimensional 
systems with chiral symmetry still remains an important 
unsolved problem. The situation can be represented by 
a single-band tight-binding model defined on bipartite 
lattices subjected to purely off-diagonal disorder like a 
random-magnetic flux with zero mean or real random 
hopping terms. The latter belongs to the chiral orthog- 
onal universality class while the former is chiral unitary. 
Due to the chiral symmetry, the model exhibits metallic 
behavior only at energy E = OJ- It is therefore of con- 
siderable interest to study the critical properties of the 
model in the vicinity of the critical point and to investi- 
gate its universality. 

The traditional finite-size-scaling analysis of disorder 
driven metal- insulator transitions, i.e., continuous quan- 
tum phase transitions at zero temperature, is based on 
two assumptions^— (i) In the vicinity of the critical point 
all variables of interest are a function of only one param- 
eter, 



v{E,L)=F{L/ti{E)), 



(1) 



where £,(E) is the energy dependent correlation length. 
Here, vicinity means that both the system size L and 
are already larger than any other typical length of 
the model and £ > L. (ii) At the critical energy E c , the 
correlation length £,(E) diverges as 



\E — E c 



(2) 



with a universal critical exponent v. Relation ^ was 
confirmed in a multitude of numerical work on disor- 
dered systems in spatial dimension 2 < d < 5 and var- 
ious physical symmetries.— As scaling variables, for ex- 
ample, the localization length; 7 - the smallest Lyapunov 
exponent^ the two-terminal conduct ance^^^ the energy 



level spacings;ii~— and the inverse participation ratioi^ 
were successfully used. 

Recently, the analysis of numerical data obtained for 
the energy dependence of the two-terminal conductance 
g on a bricklayer lattice^ which represents a generic lat- 
tice model for graphene, has led to a power-law energy de- 
pendence of the correlation length £,(E) oc where 
the critical exponent v is close to 1/3. Although this out- 
come is in agreement with previous numerical results^ - — 
for square and hexagonal lattices, it is at variance with 
the Harris criterion^! which states that v > l/d, where 
d = 2 is the cuclidian dimension of the system. More im- 
portantly, all numerical data obtained to the present date 
do not agree with theoretical predictions according to 
which the correlation length depends logarithmically on 
the energy, 

= & exp[A^ln(E /\E\)], \E\ « E , (3) 

where A is related to the longitudinal conductivity and 
Eq is assumed to be of the order of the energy band 
widthfi A possible explanation of this disagreement be- 
tween theory and numerical experiments may be that 
the energies investigated in the numerical studies, down 
to 10~ 10 so far— (in units of the hopping energy), are not 
sufficiently small in comparison to the unspecified param- 
eter Eq introduced in the theory. Thus, it could be that 
the energy interval \E\ <C Eq, where the scaling holds, 
was not reached in previous numerical studies. A second 
obstacle is the vanishing of the density of states p{E) 
which occurs at E = for hexagonal and bricklayer lat- 
tices in the presence of random-magnetic-flux disorder ;22, 
This behavior persists even in strongly disordered chi- 
ral systems so that the two-terminal conductance, which 
nevertheless turns out to be finite ~ e 2 /h at the Dirac- 
point in graphene,— is not a suitable scaling variable 
for numerical studies. Therefore, it is expedient to inves- 
tigate instead the smallest Lyapunov exponents, which 
are associated with the localization length and are not 
directly affected by the vanishing density of states. 
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In this paper, we analyze the scaling behavior of both 
the two-dimensional bricklayer and square-lattice model 
with random-magnetic-flux disorder. Using the transfer- 
matrix method for quasi-ld systems^ - — £ we calculate the 
two smallest Lyapunov exponents z\ and z 2 for energies 
very close to E = 0, not achieved in previous work. Lya- 
punov exponents are more suitable than the conductance 
since they are less sensitive to the energy dependence of 
the density of states. Another reason for using Lyapunov 
exponents is that the analysis of the conductance is far 
more time consuming, which is crucial since quadruple 
precision is necessary in our case when energies smaller 
than 10 -16 are considered. We show that in the vicinity 
of the critical point our numerical data for the Lyapunov 
exponents lead to the relation 



z 1 , 2 (E,L)=F 1 , 2 



hiL \ 



(4) 



where L is the width of the system. We prove that the 
correlation length £,(E) depends logarithmically on the 
energy (see Eq. ([3])) in agreement with the predictions 
by Fabrizio and Castclliani^ Thus, our results resolve 
the previous discrepancy between analytical theory and 
numerical calculations. 
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FIG. 1: (Color online) The smallest Lyapunov exponents 
zi(E,L) — zi(E = 0,L) (lower branches) and Z2(E,L) — 
Z2(E — 0,L) (upper branches) as a function of 1/| ln|_E|| for 
energies \E\ < 3 ■ 10~ 10 . The applied random flux strength 
is / = 0.5h/e and the width of the quasi-ld systems is in 
the range 8 < L < 64. Solid lines are quadratic fits. The 
inset shows the size dependence of z\o — z\(E = 0, L) for 
8 < L < 192. 



II. THE MODEL AND METHOD 

We study a single-band tight-binding Hamiltonian de- 
fined on a two-dimensional square lattice with nearest- 
neighbor hopping and random-flux disorder, which is in- 
troduced by complex phase factors in the transfer terms, 

I p -i8x,y-a;m,y „ \ 

T^c o xy v x _y- a J 

"h { c x,y c x+a,y + c x,y c x~a,y) , (5) 

x, y 

where c x y and c X:V denote creation and annihilation op- 
erators of a fermionic particle at site (x,y), respectively. 
For bricklayer lattices, the prime at the first sum in §5§ in- 
dicates that only the transfers along every other vertical 
bond are included. In this way, the square lattice is trans- 
formed into a bricklayer where the coordination number 
is reduced to three nearest-neighbor sites. The bricklayer 
lattice has the same topology as the honeycomb lattice 
of graphene and the Hamiltonian possesses the same 
eigenvalues zte^. The phases, which are chosen to be as- 
sociated only with the vertical bonds in the y-direction, 

&x,y;x,y+a = 6x+2a,y;x+2a,y+a — ^J^4>x,y, are defined by 

the random magnetic flux 4> x , y , which is uniformly dis- 
tributed —f/2 < 4» X: y < f/2 with zero mean and disorder 
strength < f/(h/e) < 1. The random magnetic flux 
is pointing perpendicular to the 2d lattice and periodic 
boundary conditions are applied in the y direction. In 
contrast to diagonal disorder, this random flux preserves 



the chiral symmetry for both the square and bricklayer 
lattices. We fix the units of energy and length scales by 
the nearest-neighbor hopping energy V = 1 and the lat- 
tice constant a = 1, respectively. The disorder strength is 
taken to be / = 0.5 h/eior the bricklayer and / = 1.0 h/e 
for the square lattice. 

We use the transfer matrix method^ and collect nu- 
merical data for the two smallest Lyapunov exponents 
zi(E,L) and z 2 (E,L). For the system width L and 
length L x 3> L, we calculate the transfer matrix M = 
Yif* Mi; and extract the two smallest Lyapunov expo- 
nents. The relative uncertainty e(E,L) of our data is 
2 • 10 -3 for larger widths L = 192 and L = 160, and 
decreases down to 10~ 4 for the smallest L = 8. This re- 
quires the length of the quasi-ld systems L x to be in the 
range ~ 10 8 to 10 9 . Since we expect that scaling occurs 
only in the vicinity of the E = critical point, we con- 
sider energies as small as possible, down to the point of 
\E\ = 10 -34 , at least for L < 64. This requires to perform 
the calculations with quadruple numerical precision. 

The specific symmetry of the model provides us with 
an independent test of the accuracy of our data. Due 
to the chiral symmetry, the spectrum of Lyapunov ex- 
ponents must be degenerate at the band center for all 
L, 

Zl {E = 0,L) = z 2 {E = 0,L). (6) 

Deviations from E = remove this degeneracy but the 
average value, [zi(E) + z 2 (E)]/2, equals to z\(E = 0) for 
small values of E. 
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FIG. 2: (Color online) The size dependence of the coefficients 
d and di, given by (|lip and (|12p . Solid lines are linear fits 
with slopes 2.24±0.26, 2.23±0.26 (ci and di), and 3.92±0.38, 
4.03 ± 0.40 (c 2 and d 2 ). The width is in the range 16 < L < 
192 and E = 0.8. Only data for \E\ < 3 x 10 -10 were used. 
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FIG. 3: (Color online) Finite-size scaling of the two lowest 
Lyapunov exponents for various widths L and restricted en- 
ergies \E\ < 3 • 10~ 10 . The inset shows the //-dependence 
of the parameter C(L). The solid line is the power-law fit 
confirming that C(L) oc (InL) 1934 . 



III. BRICKLAYER LATTICE 

The energy dependence of the two Lyapunov expo- 
nents is plotted in Fig. Q] for various system widths L 
of the bricklayer. Our data confirm that Z\ and z 2 are 
analytical functions of the variable \j\ In Therefore, 
we approximate their energy dependence by the Taylor 
expansion 

zi{ X , L) - Zl (0, L) = co(L) + c x (L) X + c 2 (L) x 2 , (7) 
where 



\HEo/\E\)\' W 

and by a similar expansion z 2 (x, L) — z 2 (0, L) = d 2 (0, L) + 
di(L)x + d 2 (L) X 2 for the second Lyapunov exponent. 
Comparing with Eq. (2), we conclude that ln(£ /a) <C 
A[ln(£'o/|i?|)] 1 / 2 so that in what follows we consider £ ~ 
a. The expansion coefficients Cj and di are determined 
numerically. The L-dependence of the coefficient cq and 
do shows finite size corrections. For the bricklayer, we 
found that cq and do depend only weakly on L provided 
that L > 16 (data are shown in the inset of Fig. Q}. For 
instance, we obtain that z\(E = 0, L = 16) = 1.5498 ± 
0.0003 and z ± (E = 0,L = 192) = 1.557 ±0.002. 

To estimate the energy Eq, we first minimize the ex- 
pression 

v [z(E,L)-F(E,L)} 2 

j£ [z(E,L)e(E,LW ' 1 ' 

where 

F(E, L) = (3 + A(lni) Ql X + /3 2 (ln L) a - X \ (10) 

with respect to parameters a, /?, and E . We found that 
X possesses a minimum when 2.2 < u\ < 2.5, 3.8 < a 2 < 



4.1, and 0.1 < Eq- It was not possible to obtain a better 
estimation of the critical parameters from this procedure 
since small variation of Eo can be compensated by small 
change of fi 2 and a 2 . Instead, in a more accurate analysis, 
we fit our numerical data for z\ and z 2 to the quadratic 
expansion (JT))- Fig. [2] shows the L-dependence of the 
coefficients C\ i and d\ 2 . Our data confirm the assumed 
logarithmic behavior of all coefficients occurring in the 
expansion 

ci(L), di(L) oc [ln(L)p (11) 

and 

c 2 (L), d 2 {L) oc [ln(L)p , (12) 

where a,\ and a 2 are close to the anticipated values 2 and 
4, respectively. 

Fig.[3]shows another test of the scaling of the Lyapunov 
exponents. Following the conventional scaling method^ 
we re-scaled the horizontal axis for the data shown in 
Fig. [T]by the parameter C(L): \ xC{L). The such 
obtained C(L) gives us directly the required scaling be- 
havior as shown in the Fig.[3J The data for both z\ and z 2 
scale to one universal curve. The inset to Fig. [3] confirms 
the expected power-law relation C(L) cx (lni) 2 . 

To obtain a quantitative estimation of the energy Eq, 
we repeated the scaling analysis shown in Fig. [3] for var- 
ious Eq. Although we recovered the scaling behavior 
similar to that shown in Fig. [1] (data not shown) , the in- 
dependence of the parameter C(L) depends on the choice 
of E . As shown in Fig. fH C(L) oc (lnL) K with the ex- 
ponent k decreasing when Eq increases, converging to 
k ps 2 for Eo > 1. We conclude that the energy Eo is of 
the order of unity in our bricklayer model. 

Finally, we plot in Fig. [5] the two smallest Lya- 
punov exponents as a function of a single parameter 
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FIG. 4: (Color online) The function C(L) obtained by a scal- 
ing analysis of the data for four different values of the energy 
Eo- Solid lines are power-law fits C(L) oc (\nL) K with the 
value of the exponent re given in the legend. 



(\nL) 2 /ln(E /\E\) with E = 0.8. To reduce the finite 
size corrections, we subtract from the data the values 
Z\{E = 0.L) and z 2 (E = 0,L), respectively. All data 
collapse onto a single curve. 



IV. SQUARE LATTICE 

Another possibility to check for logarithmic scaling at 
chiral quantum critical points is the numerical analysis of 
a simple square lattice^ - — We calculated Z\{E,L) and 
z 2 (E, L) for L even and Dirichlet boundary conditions in 
the transverse direction. We found that the scaling anal- 
ysis is more difficult than for the bricklayer. First, the 
finite size effects are more pronounced (see lower inset 
in Fig. [SJ. We can eliminate them, at least partially, by 
subtracting the value z(E = 0,L) from z(E,L)£2- Sec- 
ond, in the unperturbed model the van Hove singularity, 
which appears at E = compared with E = ±1 for the 
bricklayer, may spoil the scaling analysis. More impor- 
tantly, following the same procedure as for the bricklayer, 
we found that the function X given by Eq. possesses 
a minimum only for small Eq ~ 10~ 4 although the en- 
ergy band widths are about the same for both lattices. 
Since the energy E must be much smaller than Eq, we 
had to restrict our analysis to energies \E\ < 10~ 20 . For- 
tunately, in such a narrow energy interval, we can ne- 
glect the quadratic term in the Taylor expansion ([7]). As 
shown in Fig. [6l both Z\ and z 2 are linear functions of x 
when L < 32. This enables us to estimate the exponent 
a\ from the analysis of the size dependence of the slope, 
Ci(L) oc (lnL) Ql . This analysis is independent on both 
the choice of Eo and finite size effects, provided that the 
latter do not depend on the energy. 

However, the fit turns out to be rather unstable to 
small changes of the data ensemble. First, the interval 
of x is very narrow and almost all data points are ac- 
cumulated in the right part of this interval. Therefore, 
the resulting fit is very sensitive to the exact value of 




( In L ) 2 / In (E /IEI) 

FIG. 5: (Color online) The data-collapse of the first two Lya- 
punov exponents zi(E, L) — zio(L) and zi(E, L) — Z2o(L) plot- 
ted as a function of a single parameter (mZ/) 2 /ln(.Eo/|-E|) 
with Eo = 0.8. For clarity, the data for the second Lyapunov 
exponent are shifted vertically. 

Zi(E = 0). Second, although we calculated our data 
with high accuracy, Fig. [5] shows that this is still not suf- 
ficient for a perfect determination of the slope. To check 
the accuracy of a±, we tested various data ensembles and 
found that ct\ varies between 1.7 and 2.1. Nevertheless, 
our data for the square lattice are compatible with a log- 
arithmic scaling relation. 



V. CONCLUSIONS 

We analyzed the scaling behavior of the two small- 
est Lyapunov exponents z\ and z 2 in disordered two- 
dimensional chiral systems defined on a bricklayer and 
on a square lattice. We found that both z\ and z 2 follow 
a logarithmic scaling relation as considered by Sittler and 
Hinrichsen2£ with a correlation length proposed by Fab- 
rizio and Castellianii22 According to Ref . [23, the physical 
origin of logarithmic scaling is associated with multifrac- 
tality and local scaling invariance. To the best of our 
knowledge, the results presented above are the first nu- 
merical confirmation of the scaling relation Q , in which 
the scaling variable is given by the ratio of the logarithm 
of the system width L to the logarithm of the correlation 
length £, instead of the ratio L/£ as applied usually. This 
scaling is accompanied by the logarithmic energy depen- 
dence of the correlation length hi^(E) oc ^J\ii{Eo/\E\) 
valid for \E\ <C Eq. Our results also solve the contradic- 
tion between previous numerical work, which apparently 
did not reach the true scaling regime, and the Harris cri- 
terion. 

Two methods of the scaling analysis were used. Both 
confirm that the logarithmic scaling is observable only for 
very small values of the energy close to the chiral quan- 
tum critical point at E = 0. In order to resolve the Lya- 
punov exponents for energies down to \E\ = 10 -34 , the 
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FIG. 6: (Color online) The energy dependence of the two 
smallest Lyapunov exponents Zi(E) — Z\{E — 0) and Z2 — 
z%(E = 0) calculated for the square lattice. The applied 
random-flux strength is / = 1.0/i/e. Owing to the small value 
of Eq, we restricted the energy interval to \E\ < 10 -20 . Then, 
both zi and zi are linear functions of x- The first inset shows 
the L-dependence of the coefficient ci of the Taylor expan- 
sion (0. The exponents ot\ obtained are 1.92 ± 0.24 for z\ 
and 1.85 ± 0.25 for Zi. The plot of Zi(E = 0, L) vs. 1/lnL 
demonstrates the finite size effects which are much stronger 
than for the brick layer (lower inset). 



implementation of quadruple precision in the numerical 
algorithms was necessary. This probably explains why 
logarithmic scaling was not observed in previous numer- 
ical worki 1 ^— 



The question arises whether the same logarithmic scal- 
ing analysis can be performed also for the two-terminal 
conductance. At present this seems not possible with 
our available computing power. In our previous work^ 
we found g(E,L) = g ]n(E*(L)/\E\) for \E\ > E* = 
E* /const., but did not observe any energy and system 
size dependence of the ensemble averaged conductance 
g c ~ 4/ire 2 /h as long as the energy remains smaller than 
a certain value E* oc L~ 2 . This size dependent energy 
interval coincides with the recently observed depression 
in the density of states^ The observation of a tiny loga- 
rithmic energy dependence of the conductance, if present 
at all, would require a far more accurate numerical deter- 
mination of the ensemble averaged mean conductance. 

As shown analytically in Ref. [H, the logarithmic en- 
ergy dependence of the correlation length near E = 
is accompanied by a divergence of the density of 
states p(E) oc E^ 1 exp[-(4A\nE /Ey^}. Such a re- 
lation is, however, not found in recent numerical work 
on a unitary chiral lattice model, where the density 
of states decreases to zero when E — > 0<2£ Also, a 
different divergency exponent of the density of states 
p(E) oc ^~ 1 exp[-l/2(c|ln J E;/ J E;o|) 2/3 ] 

was derived an- 
alytically for the chiral orthogonal model . 28 ' 29 This dif- 
ference shows also up in the energy dependence of the 
correlation length. It would be very interesting to see, 
whether this subtle difference can also be observed in nu- 
merical studies on a bricklayer lattice with real random 
hopping disorder, which belongs to the chiral orthogonal 
symmetry class. 
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This does not eliminate the finite size effects completely, uncertainty, 
since the value of z(E — 0,L) is known only with some 



